% try many potential gearbox designs to find feasible solutions
% AOD 10.24.2011
clear
clc

switch '2' % 1 for user input design space, 2 for saved space
    case '1'
        % diametral pitches to be searched
        P = [2,2.25,2.5,3,4,6,8,10,12,16,20,24,32,40,48,64,80,96,120,150,200]; %all
        P = [2.5,3,4,6,8,10,12,16]; % narrow to subset coarse only
        P = [8,10,12,16,20,22]; % even narrower for quick calculation time, may be too narrow

        % design variable ranges
        N_1 = [16*2:22*4];
        N_2 = [15:20];
        N_3 = [16*2:22*7];
        N_4 = [15:20];
        P_12 = P;
        P_34 = P;

        % check run time before committing
        num_it=length(N_1)*length(N_2)*length(N_3)*length(N_4)*length(P_12)*length(P_34);
        estimateruntime = num_it/290000; % for AOD laptop
        text=['estimate run time: ',num2str(estimateruntime),' seconds'];
        disp(text);

        possible = zeros(num_it,6); % preallocate possible design matrix
        feasible = zeros(num_it,6); % preallocate potential design matrix

        i = 1;
        %tstart=tic; % record start time

        % established that minimum is definitely 16
        for P12=P_12; % for each standard diametral pitch on gear set 1
            for P34=P_34; % for each standard diametral pitch on gear set 1
                % note factor of 6 rounded up from sqrt(18) ideal ratio
                for N1 = N_1; % for range of teeth (GS1 gear)
                    for N2 = N_2; % for range of teeth (GS1 pinion)
                        for N3 = N_3; % for range of teeth (GS2 gear)
                            for N4 = N_4; % for range of teeth (GS1 pinion)
                                % standardize input for gearboxconstraint.m
                                x=[N1,N2,N3,N4,P12,P34];
                                possible(i,:)=x;
                                i = i+1;
                            end
                        end
                    end
                end
            end
        end

        [g]=gearboxconstraint(possible(1:i-1,:));
        ifeas = find(all(g'));

        %toc(tstart) % display total calculation time
        numberfeasiblesolutions = length(ifeas)
        feasible = possible(ifeas,:) % display feasible solutions
        save('designspace.mat','feasible')
        
        return
    case '2'
    AAA = load('designspace.mat','feasible');
    possible = AAA.feasible;
    
    num_it=length(possible);
    estimateruntime = num_it/290000; % for AOD laptop
    text=['estimate run time: ',num2str(estimateruntime),' seconds'];
    disp(text);
    
    [g]=gearboxconstraint(possible);

end

ifeas = find(all(g'));

%toc(tstart) % display total calculation time
numberfeasiblesolutions = length(ifeas)
feasible = possible(ifeas,:) % display feasible solutions
% display all gear ratios (debug)
% GR = (feasible(:,1)./(feasible(:,2)).*(feasible(:,3))./(feasible(:,4)))

% find minimum height gearbox
% from diametral pitch the diameter of the gears can be found
% gearbox height = d1+d2/2+d3/2 = n1/p12 + n2/2/p12 + n3/p34
h = feasible(:,1)./feasible(:,5) + feasible(:,2)./feasible(:,5)/2 + feasible(:,3)./feasible(:,6);
[min,ind]=min(h);
smallestgearbox = feasible(ind,:)
smallestgearboxsize = min